######################################################################################
######################################################################################
#Replication code for models presented in:

#Buntaine, Mark T. 2011. Does the Asian Development Bank Respond to Past Environmental Performance when Allocating Environmental Risky Financing? World Development 39(3): 336-350.

#Contact (as of 10/2011): mark.buntaine@duke.edu / 919.627.7066

#Current for R Version 2.12.0
######################################################################################
######################################################################################


#############################SETTING UP#############################
library(lme4) #Version 0.999375-37
library(betareg) #Version 2.2-3

#Entering data to R, file path must be changed accordingly
data <- read.csv("/Users/mbuntaine/Desktop/ADB_WorldDevelopment/Replication/WD2011_data.csv")


############Figure 1. MIXED-EFFECT PROBIT MODEL (GATE-KEEPING)########################

#Model Specification for Figure 1
fig1.model <- lmer(A.bin ~ ENV.rep + OVERALL.rep + GOVT.eff + USODA.prop + JapODA.prop + logGNI + loan.cat + island + nproj + A.two + A.three + A.four + A.five + A.six.plus + (1|country), family=binomial(link = "probit"), data=data)
summary(fig1.model)
##Note: Dropping ENV.rep from model does not change direction of other variables
##Note: Correlation between ENV and OVERALL rep is 0.230


############Figure 3. MIXED EFFECT PROBIT FOR NON-RISKY (GATE-KEEPING)################

#Model Specification for Figure 3
fig3.model <- lmer(C.bin ~ ENV.rep + OVERALL.rep + GOVT.eff + USODA.prop + JapODA.prop + logGNI + loan.cat + island + nproj + C.two + C.three + C.four + C.five + C.six.plus + (1|country), family=binomial(link = "probit"), data=data)
summary(fig3.model)


#############Figure 4. BETA REGRESSION (ALLOCATION STAGE)#############################

data.pos <- subset(data,amount.A>0)
data.pos$A.three.plus <- data.pos$A.three + data.pos$A.four + data.pos$A.five + data.pos$A.six.plus

#Model Specification for Figure 4 (top bars)
fig4.model.top <- betareg(amount.A.prop ~ ENV.rep + GOVT.eff + USODA.prop + JapODA.prop + logGNI + loan.cat + country.per + Azerbaijan + Bangladesh + China + India + Indonesia + Laos + Nepal + Pakistan + SriLanka + Vietnam + A.two + A.three.plus, data=data.pos)
summary(fig4.model.top)

#Model Specification for Figure 4 (bottom bars)
fig4.model.bottom <- betareg(amount.A.prop ~ ENV.rep + OVERALL.rep + GOVT.eff + USODA.prop + JapODA.prop + logGNI + loan.cat + country.per + Azerbaijan + Bangladesh + China + India + Indonesia + Laos + Nepal + Pakistan + SriLanka + Vietnam + A.two + A.three.plus, data=data.pos)
summary(fig4.model.bottom)


#############FUNCTIONS FOR COMPUTING BAYESIAN REPUTATION MODEL########################

##Functions for predictive distributions
#Function to compute the prior predictive of a beta-binomial
bb.prior<-function(a,b,sim.n,pred.n){
  pi.prior.samp <- rbeta(sim.n, a, b) 
  y.prior.pred <- rbinom(sim.n, pred.n, pi.prior.samp)
  return(mean(y.prior.pred))
}

#Function to compute the posterior predictive of a beta-binomial
bb.post<-function(y,n,a,b,sim.n,pred.n){
  pi.post.samp <- rbeta(sim.n, y+a, n-y+b) 
  y.post.pred <- rbinom(sim.n, pred.n, pi.post.samp)
  return(mean(y.post.pred))
}

data$country <- as.character(data$country)
proj.eval$country <- as.character(proj.eval$country)
#Note: for these functions to work, source data needs to be imported 

# Computing the ENV.rep variable from the predictive distributions
ENV.rep <- rep(0,nrow(data))
data.A<-data.frame(data,ENV.rep)
rep.window<-5 #This is length of moving window for evaluations to have impact
L<-nrow(data)
a<-3 #This is the a value for beta prior, based on observed population distribution
b<-4 #This is the b value for beta prior, based on observed population distribution
sim.n<-100000 #This is number of simulations to run for predictive distribution
pred.n<-3 #This is number of scale categories
for (i in 1:L){
  WIN<-subset(proj.eval, country==data$country[i] & yr<=data$yr[i] &
              yr>=(data$yr[i]-rep.window+1) & ENV.rating>=0)
  data$ENV.rep[i]<-ifelse(nrow(WIN)==0,
                      bb.prior(a,b,sim.n,pred.n),
                      bb.post(sum(as.numeric(WIN$ENV.rating)),
                                  nrow(WIN)*pred.n,a,b,sim.n,pred.n))
}


#Computing the OVERALL.rep variable from the predictive distributions
OVERALL.rep <- rep(0,nrow(data))
data<-data.frame(data,OVERALL.rep)
rep.window<-5 #This is length of moving window for evaluations to have impact
L<-nrow(data)
a<-2 #This is the a value for beta prior, based on observed population distribution
b<-2 #This is the b value for beta prior, based on observed population distribution
sim.n<-100000 #This is number of simulations to run for predictive distribution
pred.n<-3 #This is number of scale categories
for (i in 1:L){
  WIN<-subset(proj.eval, country==data$country[i] & yr<=data$yr[i] &
              yr>=(data$yr[i]-rep.window+1))
  data$OVERALL.rep[i]<-ifelse(nrow(WIN)==0,
                      bb.prior(a,b,sim.n,pred.n),
                      bb.post(sum(as.numeric(WIN$OVERALL.rating)),
                                  nrow(WIN)*pred.n,a,b,sim.n,pred.n))
}
